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J3 ■ Abstract 



We calculate the long-distance matrix elements for the decays of the T (rjb) 
and Xb {hb) states in lattice QCD with two flavors of light dynamical quarks. 

X" 

We relate the lattice matrix elements to their continuum counterparts through 
onedoop order in perturbation theory. In the case of the leading S"-wave 
matrix element, we compare our result with a phenomenological value that 
we extract from the experimental leptonic decay rate by using the theoretical 
expression for the decay rate, accurate through relative order a s . Whereas 
estimates of the leading 5-wave matrix element from quenched QCD are 40- 
45% lower than the phenomenological value, the two-flavor estimate of the 
same matrix element is close to the phenomenological value. Extrapolating 
to the real world of 2 + 1 light flavors, we find that this matrix element 
is approximately 6% higher than the phenomenological value, but that the 
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phenomenological value lies within our error bars. We also compute the color- 
singlet and color-octet matrix elements for P-wave decays. We find the value 
of the color-singlet matrix element for 2 + 1 flavors to be approximately 70% 
larger than the quenched value and the value of the color-octet matrix element 
for 2 + 1 flavors to be approximately 40% larger than the quenched value. 
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I. INTRODUCTION 

Bottomonium is a nonrelativistic system: the velocity v of the b and b quarks in 
the center-of-mass frame is much less than unity (v 2 ~ 0.1). Bodwin, Braaten and Lepage 
have shown that, within the framework of Nonrelativistic Quantum Chromodynamics 
(NRQCD), the smallness off allows one to expand the decay rates into light hadrons and/or 
electromagnetic decay products in powers of v . Each term in this velocity expansion can 
be expressed as a finite number of terms, each of which is a product of a long-distance 
(~ 1/Mi,v ) matrix element of a four-fermion operator between bottomonium states and a 
short-distance (~ l/Mj,) parton-level decay rate. Owing to the asymptotic freedom of QCD, 
the short- distance parton-level decay rate can be calculated perturbatively. 

The S'-wave bottomonium decay rates can be expressed, through next-to-leading 
order in v 2 , as 

T( 2s+1 S 2s+1 - X) = g 1 ( 2s+1 S 2s+1 ) 21mf\( 2s+1 S 2s+1 )/M 2 

+ F i e s+1 S 2s+1 )21mg i e s+1 S 2s+1 )/M* . (1) 

Similarly the P-wave bottomonium decay rates at lowest non-trivial order in v are given by 

TC s+1 Pj ->X) = Hi( 2s+1 P J )2Imf 1 ( 2s + 1 P J )/M* 

+ HsC s+1 Pj)2lmfsC s+1 S 2s+1 )/M 2 . (2) 

The /'s and g's are proportional to the short- distance rates for the annihilation of a bb pair 
from the indicated 2s+1 Lj state, while Qi, T\, Tii, and 7i 8 are the long-distance matrix 
elements.^ The subscripts 1 and 8 indicate that the bb pair is in a relative color-singlet or 
color-octet state. If one works to leading order in v in the NRQCD Lagrangian, then the 
matrix elements of the spin-singlet and spin-triplet states are equal. 



x Our quantities Hi and Hs are related to the quantities H\ and H% in Ref. by Hi = M^H\ 
and H 8 = M%H 8 . 



In earlier papers 0, we reported lattice NRQCD calculations of Qi, T\, TCi, and 
Hs for the T (r] b ) and Xb (h b ) states that made use of quenched gauge-field configurations 
with inverse lattice spacings a -1 ~ 2.4 GeV and a -1 « 1.37 GeV. We found that the value 
of the best-measured matrix element Q\ is 40-45% below a phenomenological value that 
we extracted from the leptonic width of the T and the theoretical expression for the width, 
accurate through relative order a s .f] The NRQCD Collaboration fUH had noted that at least 
part of the discrepancy is likely due to the use of the quenched approximation. The reason 
that the quenched approximation underestimates the matrix element is that the distance 
scale associated with the bottomonium bound state [order l/(M&i>)] is considerably larger 
than the scale at which the matrix elements sample the wave function (order a, which is order 
1/Mj). If we fix the lattice QCD coupling at \j(M b v) to a value that yields good agreement 
with the bottomonium spectrum, then, in the quenched approximation, the coupling at a 
will be weaker than it should be. Hence, the wave function at the origin will be too small, 
leading to a prediction for the bottomonium decay rate that is too small. 

In this paper we present calculations of the decay matrix elements for the T (j]b) 
and Xb (hi,) states that make use of gauge configurations containing the effects of 2 flavors 
of light dynamical (staggered) quarks. These calculations confirm that most, if not all, 
of the discrepancy in the previous calculations of the matrix elements was, in fact, due to 
quenching. Our results, when extrapolated to three light flavors, lead to a slight overestimate 
of the T decay rate. 

The remainder of this paper is organized as follows. In Sec. II we define the required 
matrix elements in the continuum and on the lattice and describe the lattice implementation 
of NRQCD that we use in our calculations. Section III contains an outline of the perturbative 
calculation that we use to relate the lattice matrix elements to their continuum counterparts. 



2 The phenomenological value that we quote in the present paper is based on a slightly different 
value for a, than was used in Ref. 01. 



We present our results in Sec. IV, and Sec. V contains our conclusions. 

II. MATRIX ELEMENTS AND LATTICE NRQCD 

In the leading non-trivial order in v, the NRQCD Lagrangian for the bottom quark 
and antiquark is 

£fl = * t ( A -S)* +x '( A + S x ' (3) 

where tp is the quark annihilation operator and \ is the antiquark creation operator. D t and 
D are the gauge-covariant temporal and spatial derivatives. Note that, although one can 
obtain the correct leading-order spectroscopy in the Coulomb gauge by replacing D with 
the simple (non-covariant) gradient operator, the covariant operator is needed to calculate 
the octet P-wave decay matrix element, even at lowest non-trivial order. Not surprisingly, 
Eq. (H) is just the Euclidean-time Schrodinger Lagrangian for the bottom quarks and anti- 
quarks. 

We work to leading order in v in the Lagrangian. As we have mentioned, at this order, 
the matrix elements of the spin-triplet and spin-singlet states are identical. Therefore, we 
approximate all of the long-distance matrix elements in Eqs. ([!]) and (0) as spin-singlet 
matrix elements. Using the leading-order Lagrangian, we are able to compute the order-v 2 
.S-wave matrix element J-'i, with an error of order v 4 . Note, however, that, in order to obtain 
a full relative-order-f 2 computation of the ^-wave decay rate, we would need to compute Q\ 
through relative-order v 2 . This would require relative-order-f 2 terms in the Lagrangian, in 
which case the spin-singlet and spin-triplet states would be distinguished. 

In terms of the fields x an d if}, the spin-singlet matrix elements that we compute are 

ft = (^ol^XxWSo) , (4a) 

Fi^CSol^XxH-inM'So), (4b) 

Hx = ^PiWii/l) D X ■ X j m) D ip\ x Pi) , (4c) 



H 8 =( 1 Pi\^T a xX i T a ^\ l P 1 ), (4d) 

where x 1 D i> = X^Dip — (D%) V- 

The vacuum-saturation approximation is valid for the color-singlet matrix elements 
and is accurate up to errors of relative order v 4 (Ref. |l|). In that approximation, Eqs. 



and (HcT) become 



ging? a = ( 1 So\ii> i x\o){o\xty\ 1 So), M 

Tx « ?T = ( 1 5o|^ t x|0)(0|x t (-f D) 2 Vr^o) , (5b) 

Hi « Wy 8 = CPil^Ci/Z) D x|0> ■ (0|x f (V2) D ^\ X P X ) • (5c) 

One can express vacuum-saturation values of the color-singlet matrix elements as ^^ s = 
t|^|-Rs(0)| 2 and 1~C[ S = ^r|-R' P (0)| 2 , where Rs(Q) is the radial wave function of the S'-wave 
state at the origin and R' P (0) is the derivative of the radial P-wave wave function at the 
origin |l[. These are the quantities that appear in decay rates in the color-singlet model. 
In contrast, the term proportional to 7i§ is absent in decay rates in the color-singlet model. 
7^8 is the probability of finding a bbg component in P-wave bottomonium, with the bb in a 
color-octet state. 

In our lattice calculation of these matrix elements, we transform our gauge field 
configurations to the Coulomb gauge. For this gauge choice, we can replace the covariant D 
with the non-covariant V in Eq. @. Corrections to this replacement are suppressed by v 2 . 

We employ various discretizations of the derivative operator. For the operator 7i\, 
we replace the covariant derivative D with the non-covariant finite difference 5, which is 
defined by 

Si^(x) = ~[ip(x + i) - tjj(x - i)] , (6) 

where i is the unit vector in the ith spatial direction. For T\ we employ four different 
discretizations of D 2 : 
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A (2) (non)^(x) = YH^( X + { ) + ^( x ~ [ ) ~ 2 ^( x )} , ( 7a ) 

i 

A (2) (cov)^(x) = J2 { — [Ui(x)ip(x + i) + UJ{x - i)ij)(x - i)] - 2V»(x)]| , (7b) 

^ A ( 2 )(non 2 ) X = - EKW^x] , (7c) 

i 

^ t A( 2 )(cov 2 ) X = -DWMx] , (7d) 

i 

where the covariant finite difference d is defined by 

di^(x) = ^-[Ui{x)^{x + i) - [//(x - i)V(z - i)] , (8) 

and m is the tadpole contribution to U. We adopt the definition u Q = (l^/piaq) 1 ^ 4 - 

On the lattice, we obtain such matrix elements by measuring the expectation value in 
the gluon background of a product of three operators: a source for a bb pair with the appro- 
priate quantum numbers at a (Euclidean) time — T, the appropriate four-fermion operator 
at time zero, and a sink for the bb pair at time T". For convenience, and in order to reduce 
noise, we divide this expectation value by the product of two other expectation values. One 
is the expectation value of the product of the numerator source for the bb pair at time — T 
and a point sink that annihilates the bb pair at time zero; the other is the expectation value 
of the product of a point source that creates a bb pair at time zero and the numerator sink, 
which annihilates the bb pair at time T'. This ratio is illustrated in Fig. 0. 





FIG. 1. Lattice calculation of a matrix element of a four-fermion operator. The large discs 
represent the sources and sinks; the smaller discs represent the four-fermion and point source 
operators. The lines are the nonrelativistic quark propagators. 



In the cases of Q\ and 7ii, this ratio approaches the ratio of the matrix element to 
its vacuum-saturation approximation in the limit T,T' — > oo. Hence, it gives an indication 
of the accuracy of the vacuum-saturation approximation. In the case of Tig, this ratio yields 
'Hg/'hO[ s in the limit T,T' — > oo. We obtain values for Q\ and 7i\ in the vacuum-saturation 
approximation from the relations 

E(0|x f (x, T)i/,(z, T)^(0, 0) X (0, 0)|0) ^ grexp(-E s T) (9) 

X 

and 

£<0|xt(a;, T)(-| D)^(x,T) •^ t (0,0)(-f D)x(0,0)|0) T -^? Hj s exp(-E P T) , (10) 

which follow from the fact that only the lowest-lying intermediate state with the correct 
quantum numbers contributes to the amplitude in the limit T —>■ oo. Note that we can write 

0r = a 2 p , (11) 

where one factor of a p is from the point source and the other is from the point sink. If we 
replace the point source by another (extended) source, the coefficient of the exponential is 
of the form a p a x , while if we use this new extended-source operator for both source and sink, 
the coefficient is a x . Thus, introducing an extended source which has a greater overlap with 
the ground state gives us an alternative method of extracting a p and, hence, Q\. Similar 
comments hold for 7i\ and T\. We calculate the .TVs from 

-E x (0\x(x,TyAW(*)ij(x,T)S(0)\0) t^oo .^T 



E,(0|x(:r,T)t^(x,T)S(0)|0> ^ 



VS ' 



(12) 



where A( 2 )(*) denotes any of the discretizations of D 2 in Eq. (|7|), and 5(0) is any source 
with a finite overlap with the lowest S"-wave state on time slice 0. 

In order to evaluate these matrix elements, we must calculate bottom-quark propaga- 
tors G(x; y) on the lattice. Following Lepage et al. 0, we calculate the retarded propagator 
G>(x, t; 0) by iterating the equation 
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G r (x,x + 1; 0) = (1 - H /2n) n Ul >X0 (l - H /2n) n G r {x,x ; 0) + 5 x , 4 +i,o , (13) 

setting G(x,a;o;0) = for x < 0. In Eq. flU), H = -A^/2M - h , A® is the gauge- 
covariant discrete Laplacian, which is given by the expression in Eq. (|7b|) with uq set to unity, 
ho = 3(1 — Mo) /Mo, and M$ is the bare bottom-quark mass. We note that our bare bottom- 
quark mass is defined to be u times that of Ref. 0. The value two for the discretization 
parameter n turns out to be adequate for our calculations. 

An expression that is similar to Eq. ( pll) exists for the advanced propagator G a . The 
relation G r (x;y) = G' a (y;x) makes it possible to rewrite amplitudes, interchanging sources 
and sinks. Such a rewriting allows one to start all propagator calculations from a noisy 
(point or extended) source, rather than a point source and, thereby, to reduce both the 
statistical error and the number of calculational steps. 

III. THE RELATIONSHIP BETWEEN LATTICE AND CONTINUUM MATRIX 

ELEMENTS 

We wish to relate our lattice results to the continuum [modified minimal subtraction 
scheme (MS)] matrix elements that are used in phenomenology. Lattice matrix elements and 
continuum matrix elements differ only in the choice of ultraviolet regulator. Furthermore, 
a change of ultraviolet regulator is dependent only on the large-momentum (short-distance) 
parts of an amplitude. Consequently, asymptotic freedom allows us to compute the short- 
distance coefficients that relate the lattice matrix elements to the continuum matrix elements 
in a perturbation series in the strong coupling a s . The short-distance coefficients are in- 
dependent of the hadronic state. Therefore, for purposes of computing the short- distance 
coefficients, we choose, for convenience, to evaluate the operators in free QQ states. 

We can expand the lattice-regulated matrix element of an operator in terms of 
continuum-regulated matrix elements of a complete set of operators: 



(O i ) L = J2c ij (O j ) c , (14) 

3 

where the C{j are the short- distance coefficients, (O) is the matrix element of the operator O 
in a free QQ state, and the subscripts L and C indicate the lattice- and continuum-regulated 
matrix elements, respectively. The matrix elements and short-distance coefficients can be 
expanded in perturbation series: 

(O t ) L =(O l )^+a s (O l )i ) + ---, (15a) 

(O l ) c =(O l )P+a s (O l )£ ) + ---, (15b) 

Cij= c\f + a s c\f H . (15c) 

For simplicity, we use the same definition of a s and the same scale for a s in all three 
expansions in Eq. (|TJ). 

At zeroth order in the perturbation series, the momentum-space expression for a 
lattice operator is equal to the momentum-space expression for the corresponding continuum 
operator, plus terms of higher order in the lattice spacing a times the momenta. Therefore, 

4? = 1, (16) 



and 



4 0) = for Dim O, < Dim { , (17) 

where Dim O is the mass dimension (or, equivalently, order in v) of the operator O. For the 
operators that we consider in this paper, 

i 0) = 0for^j. (18) 

Since our lattice NRQCD action is accurate only to leading order in v, only the 
following mixings can be treated consistently: Qu, into Q\c an d T\c, J~il into T\c-, Tin 
into Tt\c and TCsc, and TCsl into Hgc and TCic- Therefore, we need consider, at most, two 
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operators in the expansion ([14]). Then, using Eqs. ([15|), flT6[), and (|T8|), we equate the terms 
of order aj in Eq. ([L4]) to obtain 

(Oi) 1 ? ~ (Ot)® = <S\0^ + cffiO^P for i + j, (19) 

where no sum over j is implied. The quantities on the left side of Eq. flT9p are computed 
in perturbation theory. We determine the short- distance coefficients on the right side of 
Eq. ([l^) by expanding the quantity on the left side of Eq. (|19"D in powers of the external 
QQ 3-momenta and by choosing free QQ states with particular color (and, in general, spin) 
quantum numbers. 

In the expansion of the quantity on the left side of Eq. ( |I9"D in powers of the external 
QQ 3-momenta, the various terms are infrared finite, to the extent that the behavior of 
the integrand in the lattice matrix element matches the behavior of the integrand in the 
continuum matrix element at small loop momentum. The expansions for the various mixings 
that we have mentioned above yield, at most, a linear infrared divergence in the lattice 
and continuum matrix elements. Since our lattice action (and, implicitly, our continuum 
action) are accurate to leading order in v 2 , those divergences cancel between the lattice and 
continuum matrix elements on the left side of Eq. fll9|) . 

In general, infrared divergences in differences between lattice and continuum matrix 
elements cancel, provided that one works consistently to a given order in v. This means that, 
in order to compute coefficient of the mixing of a lattice matrix element into a continuum 
matrix element of relative order v n , one must employ lattice and continuum actions that are 
accurate to relative order v n . Then, the small-loop-momentum behaviors of the lattice and 
continuum contributions on the left side of Eq. ([19]) will be the same, and infrared divergences 
in the mixing coefficient will cancel. On the other hand, one should not compute operator 
mixings that exceed the accuracy in v of the action. For example, since we use actions of 
leading order in v, we do not compute the one-loop correction to the mixing of Qn, into T\c 
(relative order v 2 ). If one were to carry out such a computation, using the leading-order 
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NRQCD lattice and continuum actions, then the expression for the one-loop correction would 
contain a cubic leading infrared divergence in both the lattice and continuum contributions 
on the left side of Eq. ([19|). The cubic leading divergence would cancel, but, owing to the 
absence of order- 1> 2 terms in the action, a linear sub leading divergence would persist. 

Applying Eqs. (|T^D, flTBD, and flT8| ) to the operator matrix elements that we consider 
in this paper, we obtain 

Gil = (1 + e)Gi , (20a) 

J- 1L = (l+ 7 )^i + 0^i, (20b) 

and 

H 1L = (1 + t)Hi + kH 8 , (21a) 

H 8L = (I + r])H 8 + (Hu (21b) 

where we have dropped the subscript C on the continuum matrix elements, and the coeffi- 
cients e, 7, 0, l, k, r\ and ( are of order a s . It turns out, in an explicit calculation, that the 
coefficient k actually vanishes in order a s . Details of the calculations of these coefficients 
will be given elsewhere ||. 

We note that the perturbation series for the MS continuum short- distance coefficients 
that relate matrix elements to physical quantities contain renormalon ambiguities. The MS 
continuum operator matrix elements contain compensating ambiguities, and, so, the phys- 
ical quantities are ambiguity free ||. In contrast, the lattice operator matrix elements and 
the short- distance coefficients that relate them to physical quantities are free of renormalon 
ambiguities 0. Consequently, the perturbation series that relate the lattice and the MS 
continuum operator matrix elements contain renormalon ambiguities. At the one-loop or- 
der to which we work, the factorial growth of the series associated with the presence of 
renormalons is unimportant. However, because the series that relate the lattice and the 
MS continuum operator matrix elements (and the series that relate physical quantities and 
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the MS continuum operator matrix elements) ultimately fail to converge, the value of an 
MS continuum matrix element is meaningful only if one specifies the order in perturbation 
theory that is employed in computing it. 

IV. RESULTS 

A. Lattice Computation of the Matrix Elements 

For the lattice calculations, we use gauge configurations generated by the HEMCGC 
Collaboration JTIJ with two flavors of light dynamical staggered quarks on a 16 3 x 32 lattice 



at /3 = Q/g 2 = 5.6. We use all 399 configurations with light-quark mass m = 0.01 (in 
lattice units) and 200 configurations with quark mass m = 0.025. As we have already 
mentioned, we follow Lepage et al. in choosing u$ = [|Tr[/pi aq ] 1//4 as our definition of 
the tadpole contribution to U. Our measurements yield u = 0.866985(11) at m = 0.01 
and Mo = 0.866773(12) at m = 0.025. Since these are so close, we use 0.866859 for our 
perturbative calculations. We choose our bare bottom-quark mass to be M b = 1.56 ~ 1.80%, 
where 1.80 is the value chosen by the NRQCD Collaboration [^] to yield the best fit to the 
T-Xb and T-T' mass splittings. 

To calculate the required matrix elements, we first gauge fix our configurations to the 
Coulomb gauge. We then generate the advanced and retarded bottom-quark propagators 
from a stochastic estimator to an S"-wave point source, a stochastic estimator to an S"-wave 
Gaussian source, and a stochastic estimator to a P-wave point source for each color on 
each time slice. The width of the Gaussian source is chosen to be 2.5 in lattice units, 
which is approximately the radius of T or ^. From these we calculate the S- and P-wave 
bottomonium propagators, with both point and Gaussian sources and sinks, and the matrix 
elements of Fig. [I] and Eq. (|12|)- Because the extended source has a larger overlap with the 
ground state than does the point source, we extract a 2 , from fits of the propagator with an 
extended (Gaussian) source and sink to the form a^exp(— ET) for large T. We then calculate 
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the ratio a p /a x from fits of the ratio of the propagator with extended source and point sink to 
the propagator with extended source and sink. Finally, we extract Qi(Ti.i) = 2<2q. [An extra 
factor of two appears here relative to Eq. ( [TT]) because, owing to the spin independence of the 
lattice action at leading order in v, we compute propagators for only a single spin component.] 
In the case of Qi, the direct extraction from the point-source-point-sink propagator gives 
a result that is consistent with this indirect method. However, for 7i\ the point-point 
propagator is very noisy and shows no sign of a plateau in the effective wave-function plot. 
In this case the indirect method is required. Figure § shows the effective wave function as 
a function of T for the S-wave extended-extended propagator. Figure |3] shows the ratio of 
the S'-wave extended-point propagator to the S'-wave extended-extended propagator as a 
function of T. 

Our estimates of J-"i S /Qi S from the various discretizations of D 2 are obtained from 
fits to the propagator ratios of Eq. ( |I2| ) for the extended source. Our point-source results 
for these ratios are completely consistent with the extended-source results. 

Finally, we extract the ratio 7is/T-Cf s from the quantity represented in Fig. [I], where 
the 4-point vertex denotes the octet operator of Eq. (13). Its value for the case of a point 



source and sink is plotted in Fig. [|. We consider fits over the ranges T\ < T, T" < T 2 , for all 
choices of 7\ and T 2 , excluding overlaps. From these we choose a "best" fit, i.e., one with 
a good confidence level, small error, and a reasonably large range T 2 —T\. The chosen best 
fit is over the range 2-12 and has a confidence level of 40%. It yields a value 'Hsl/'Hil — 
0.01565(8). In comparison, the fit with the highest confidence level (99.8%) is over the 
range 6-8 and yields a value TC%l/'Hil = 0.01540(16), which is in agreement the selected 
fit. The results for the extended source are consistent the results for a point source, but the 
plateau occurs roughly one unit later in T, T', and the "data" are noisier. We estimate the 
systematic error in Hsl/'Hil by examining the entire plateau, both for the point-source data 
and for the extended-source data, and determining the range of fluctuations in the region in 
which the signal-to- noise ratio is appreciable. 
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FIG. 2. Effective a\ as a function of T for S'-wave bottomonium. 

Note that there is clear evidence for a plateau in 'H&l/'H-il f° r T,T' > 1, and, so, 
we are justified in assuming that the asymptotic behavior occurs for relatively small T, T", 
where the signal-to-noise ratio is relatively good. An analysis of the effective wave function 
for the P-wave state shows a plateau that starts at T m 10. However, we observe in the 
ratios that we use to calculate T\, in which we have decent signals out to T = 31, that the 
plateau can start much before effective masses and effective wave functions indicate that 
one has obtained a pure state. This is also the case for the data that we present later on 
G\/GY S - Presumably, the early onset of a plateau in these ratios indicates that their values 
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S— wave "wave functions" point/extended 
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FIG. 3. Ratio of the extended-point propagator to the extended-extended propagator as a 
function of T for S-wave bottomonium. 

are not very different for the IP and 2P states (and IS" and 2S states). 

The results for these lattice matrix elements are given in Table |. The first error 
bar is statistical. The second error bar is a combination of our estimate of the systematic 
error that arises from our choice of fits and our estimate of the uncertainty that arises 
from the fact that the propagators have not reached their asymptotic forms in the region 
of measurement. We note that the dependence on the light-quark mass is weak. For this 
reason and for the reason that we have fewer configurations at the higher light-quark mass, 
we have not calculated the P-wave matrix elements at m — 0.025. 
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FIG. 4. TIsl/T~(-il as a function of T and T" for point sources/sinks. Error bars have been 
suppressed to make the graph more easily readable. The suffix p indicates that we used a point 
source and point sink. 

B. Lattice-to-Continuum Conversion 



First, let us present the one- loop results for the coefficients that relate the lattice 
matrix elements to their continuum (MS) counterparts. These coefficients were defined 
and the method for their calculation was outlined in the Sec. [III]. The loop integrals were 
evaluated numerically, using the adaptive Monte Carlo routine VEGAS |iT|. The values of 
the coefficients, in lattice units (a = 1), are presented in Table [TI]. These values depend on 
the value of the bottom-quark mass in lattice units. However, as we have already discussed, 
we take the bottom-quark mass, in lattice units, to be the same at both the values of the 



17 



m=0.01 



m=0.025 



Gil 
FiL(non)/gi L 

J 7 il(cov)/Q 1l 
^ii(non 2 )/^iL 

Ti L {co\ 2 )/Qil 

7~t.lL 

7~(-8l/7~(-il 



0.20479 ± 0.00036 ± 0.0028 
1.53074 ± 0.00049 ± 0.001 

0.99667 ± 0.00037 ± 0.0005 

0.86310 ± 0.00027 ± 0.0005 
1.23209 ±0.00022 ±0.001 

0.02016 ± 0.00078 ± 0.0014 
0.01565 ±0.00008 ±0.001 



0.21265 ±0.00061 ±0.0028 
1.54816 ± 0.00069 ± 0.0005 
1.00740 ± 0.00049 ± 0.0005 
0.87226 ± 0.00037 ± 0.0003 
1.23961 ± 0.00029 ± 0.0005 



TABLE I. Lattice bottomonium decay matrix elements for light quark masses m = 0.01 and 
m = 0.025. 

light-quark mass m that we use. Then, with the exception of £, the coefficients in Table |T.( 
depend on the light-quark mass only through the scale of a s , which is proportional to a, 
since a depends (weakly) on m. ( has additional dependence on a and, hence, on m, since 
it contains a term that is proportional to ln(/ia), where /x is the NRQCD factorization scale. 
We take /i = 4.3 GeV, which is close to M 6 (MS). 

We convert the lattice matrix elements to continuum matrix elements using the 
formulae of Eqs. (|20|) and (j2ip . Here, we choose a s = ap(l/a), where ap is defined 
in Ref. |7]]. To convert to physical units, we use a -1 = 2.44 GeV for m = 0.01 and 



a -1 = 2.28 GeV for m = 0.025, as determined by the NRQCD Collaboration |7| from the 
T~Xb mass splitting. Then, the required values of a s are ap(2.44 GeV) = 0.2941 ± 0.0070 
and «p(2.28 Gev) = 0.3056 ±0.0076, respectively. Our continuum matrix elements are given 
in Table |TJ. The first two error bars arise from the statistical and systematic uncertainties 
in the lattice calculation. The third error bar is our estimate of the uncertainty from un- 
calculated two-loop corrections to the coefficients in Table [nj This uncertainty is estimated 
as the greater of a s times the one-loop contribution and a 2 s times the tree-level coefficient. 
Clearly this is the dominant uncertainty. In the case of J-\, this uncertainty is magnified 



Coefficient 



Value 



e 

7(11011) 

7(cov) 

7(110112) 

0(non) 

(f)(cov) 

</>(non 2 ) 

1 

V 

C(m = 0.01) 

C(m = 0.025) 



-0.4387 a s 

-0.9622 a s 

-2.543 a s 

-0.9489 a s 

4.822 a s 

3.729 a s 

3.078 a s 

-1.232 a s 

0.05484 a s 

-0.01285 a s 

-0.01680 a. 



TABLE II. Coefficients relating lattice and continuum matrix elements. The different versions 
of 7 and 4> relate to the different discretizations of D 2 . 



because the left side of Eq. fl20b|) is very close in size to the second term on the right side 
of Eq. ( |20b|) . Consequently, our calculation of T\ is very imprecise. The lattice operator 
matrix elements JF^cov), ^(non), and jF 1 (non 2 ) all yield values of T\ that are consistent 
with zero. Furthermore, the error bars in each case are larger than the differences between 
the central values. The operator matrix element .7-1(110112) yields the smallest uncertainties, 
and it is the value that derives from this matrix element that we report in Table |T|. 



There are some additional uncertainties that are not included in Table III. One is the 



uncertainty that arises from the uncertainty in the NRQCD Collaboration's determination 
of a [0]. We include only the statistical uncertainty in a -1 in our calculation. The NRQCD 
collaboration also reports an order-a 2 uncertainty and an order-f 4 uncertainty. The former 
is equivalent, in NRQCD, to the order-u 2 uncertainty, which we estimate later. The latter we 
ignore in comparison with the order-f 2 uncertainty in our calculation. The uncertainty in a 
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translates into an uncertainty in Q\ of approximately 8% at m = 0.01 and 19% at m = 0.025. 
In the case of Tii, it leads to an uncertainty of approximately 14% at m — 0.01. Another 
uncertainty arises from the neglect of corrections of higher order in v 2 in the action. These 
are nominally of order 10%, but are expected to be much smaller for spin-averaged quantities. 
(Of course, in order to obtain a spin-averaged value of Qi, one would need to observe the rjb 
and measure its decay width into 77.) Finally, there is the effect of order-a 2 corrections in 
the gluon and light-quark sectors, which cause appreciable flavor-symmetry breaking at the 
values of the lattice spacing that we use. These effects could best be estimated by repeating 
the calculation at a different value of (3. 

C. Phenomenological Value of the Matrix Element 

We obtain a phenomenological estimate for Qi from the leptonic decay width of T 



m 



Here, we use M 6 (pole) = 5.0 ± 0.2 GeV |J, a ( M b) = 1/132, a s (M b ) = 0.212, and 
r(T — > e + e~) = 1.32 ± 0.05 Gev [14]. The value of Q\ given in Table |l| includes only 



the experimental uncertainty. 

In extracting the phenomenological value of Q\, we have not included the relative- 
order-a 2 correction to T(T — > e + e _ ) fl5| . It would be inconsistent to include this correction 
without also including the order-a 2 corrections to the short- distance coefficients that relate 
the lattice operator matrix elements to the continuum ones. The relative-order-a 2 correction 
to T(T — > e + e~) contains a large dependence on the NRQCD factorization scale fi. If we 
did include this correction in our extraction, then the phenomenological value of Q\ would 
range from 3.76 GeV to 8.77 GeV as the [i ranges from 1 GeV to Mb. This large /1 
dependence and the large size of the correction at // = Mb would seem to indicate that 
the uncertainty in the phenomenological value may be close to 100%. However, experience 
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with one-loop corrections to quarkonium decay processes suggests that such large corrections 
may be canceled by the order-a 2 corrections to the short- distance coefficients that relate the 
lattice operator matrix elements to the continuum ones. Certainly, the large /i dependence 
in the relative-order-a 2 correction to L(Y —>■ e + e~) would be compensated exactly by a large 
/x dependence in the order-a 2 corrections to the short- distance coefficients. 

The uncertainty in the phenomenological value of Q\ that arises from the uncertainty 
in the value of M^ is about 8%. This is negligible in comparison with the uncertainty 
associated with the perturbation expansion. Given present theoretical uncertainties, it is 
not yet possible to extract T\ from experiment. 

Tii and TCg are related to the Xb decay widths, which have not yet been determined 



in experiments. Large corrections to the perturbation series JT1J] are likely to be important 
sources of uncertainty in the determination of these quantities, once experimental data 
become available. 



Calculation (nj = 2) 
Lattice Units Physical Units 



Phenomenology 



m = 0.01 






Qi 


0.2351(4)(32)(240) 


3.416(6) (47) (340) GeV 3 


TxlQx 


-0.8 — 0.3 


-5 — 2 GeV 2 


m 


0.032(1) (2) (5) 


2.7(1)(2)(5) GeV 5 


Hs/Hi 


0.01354(5)(63)(390) 


0.002275(9) (105) (660) GeV" 2 


m = 0.025 






Qi 


0.2456(7)(32)(270) 


2.911(8)(38)(320) GeV 3 


TxlQx 


-0.9 — 0.3 


-4.7 — 1.5 GeV 2 



3.86(14) GeV 3 



3.86(14) GeV 3 



TABLE III. Continuum MS bottomonium decay matrix elements from our lattice calculations 
with two dynamical light quarks (rif = 2) and, for comparison, a phenomenological value of Q\. 
The error bar on the phenomenological value of Q\ does not include the theoretical uncertainty. 
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D. Extrapolation to Physical Light-Quark Values 

We use linear extrapolation methods to estimate the calculated matrix elements at the 
physical values of the light-quark masses and at the physical number of light-quark flavors. 
Extrapolating to m = 0, we find that Q\ = 3.75(1) (8) (38) GeV 3 . To extrapolate to the 
physical values of the light-quark masses, we use the HEMCGC light-hadron spectroscopy 
measurements on the gauge configurations that we employ QTOJ to estimate that one-third 
the mass of the strange quark is approximately 0.0071, in lattice units. Then, we extrapolate 
Q\ to this value of m. Note that, since we are using linear extrapolations in both m and 
the number of flavors, this procedure yields the same result as would setting m s = 0.02 and 
m u = Md = 0. Finally, we use our results for quenched QCD at (3 = 6.0 || to extrapolate to 
three light-quark flavors, obtaining Q\ = 4. 10(1) (9) (41) GeV 3 . This result is approximately 
6% higher than the phenomenological value. Similarly, extrapolations to three light-quark 
flavors (with no extrapolation in m) yield Ji\ ~ 3.3 GeV 5 and TisfHi ~ 0.0018 GeV -2 . 

E. Tests of the Vacuum- Saturation Approximation 

Our lattice calculations permit us to test the validity of the vacuum-saturation ap- 
proximation for Qi and Tii. NRQCD predicts that 

Gi/Gr = 1 + 0(v 4 ) (23) 

n 1 /nr = i + o{v 4 ). (24) 

Note that, although our lattice action is accurate only to leading order in v, it does contain 
interactions of relative order v 2 , which arise through the terms proportional to the gauge 
field in the covariant derivatives. These interactions allow for the spin-independent emis- 
sion of transverse gluons, which produces the leading correction to the vacuum-saturation 
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approximation.^ 

In order to test the relations (^3) and (plfi, we need to observe a plateau in ratios of 
the form of Fig. |1[ We have measured these ratios with both point and extended sources. 
The detailed method of analysis is similar to that for the ratio / H%l/7~(-il described above. 
Because Q\ oc Q\l and Tii oc Hil, up to corrections of relative order a s , we use the lattice 
quantities to evaluate the ratios in Eqs. (EBf) and (P|). We find that 



gjg™ = i.ooi7(i) (25) 

T-d/H^ = 1.0049(2) . (26) 

These results are consistent with v 2 being of order 0.1 and justify our use of the vacuum- 
saturation approximation in computing matrix elements. Figure |5] shows the plateau in the 
ratio of lattice matrix elements Gil/Gil- We note that the plateau is reached for T,T' > 1. 



F. The Nonrelativistic Energy 

For comparison with the work of the NRQCD Collaboration, we give our estimate 
for their "energy" E NR || , which is related to the E$ of Eq. (||) by 

E NR = E s + 2\nu . (27) 

For m = 0.01 we obtain E NR = 0.4841(2), and for m = 0.025 we obtain E NR = 0.4901(3) 
(statistical errors only). This is to be compared with the value E^ R = 0.493(1) (m unspec- 



3 The inclusion of covariant derivatives does not, of course, generate all of the relative-order-t; 2 
contributions to the matrix elements we consider. In particular, the action that we use contains 
none of the spin-dependent terms that distinguish states with the same orbital angular momentum 
but with different total angular momentum (such as the T and the %). 
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FIG. 5. The ratio Gil/GYl as a function of T and T", where the subscript x indicates the use 
of an extended source. The error bars have been suppressed for clarity of presentation. 



ified) that was obtained by the NRQCD Collaboration ||TB| using an action that is accurate 
through next-to-leading order in v 2 . 



V. DISCUSSION AND CONCLUSIONS 

We have measured, in lattice simulations with two light-quark flavors, the matrix 
elements of leading order and next-to-leading order in v 2 that mediate the decays of the T 
and the color-singlet and color-octet matrix elements of leading order in v 2 that mediate 
the decays of the Xb states. We have also computed the relations between the lattice matrix 
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elements and continuum (MS) matrix elements to one-loop accuracy and have converted our 
lattice results to continuum results. We have extrapolated these results to the case of three 
light-quark flavors. 

Our measured value of the lowest-order T-decay matrix element, extrapolated to 
three flavors, is Q\ = 4.10(1) (9) (41) GeV 3 . This is in good agreement with the phenomeno- 
logical value Q\ = 3.86(14), which we extracted from the experimental value for the leptonic 
decay width by using the perturbative-QCD expression for the leptonic decay width, accu- 
rate through relative order a s . This contrasts with the value for Q that we obtained in the 
quenched approximation ||, which is 40-45% lower than the phenomenological value. The 
large size and large scale dependence of the order-a;? correction to the leptonic width fj~5 



suggest that the theoretical uncertainty in the phenomenological value of the matrix element 
may be quite large. However, one can utilize the order-a^ correction consistently only when 
one has incorporated two-loop corrections into the coefficients for the lattice-to-continuum 
conversion. Then the scale dependence would be compensated exactly, and some of the large 
correction would likely be canceled. 

At present, the most that we can say about the next-to-leading-order T-decay ma- 
trix element T\ is that the continuum (MS) matrix element is probably negative. Notice 
that, although the lattice matrix element is strictly positive, the subtractions involved in 
converting it to a continuum matrix element can, and apparently do, change its sign. If this 
result is maintained at higher orders in perturbation theory, then it is clear that any simple 
potential model, which must, of necessity, give a positive result, cannot yield the correct 
values for such higher-order matrix elements. 

We find that the value of the color-singlet P-wave matrix element 7i\ for three light- 
quark flavors is roughly 70% higher than the quenched value. The color-octet contribution 
to the P-wave decays, which is mediated by 7ig, arises from a distinctive QCD effect: the 
process in which the bb color-singlet P-wave state fluctuates into a bb color-octet S"-wave 
state plus a soft gluon. Such a contribution is absent in simple potential models. We find 
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that the value of 7i% for three light-quark flavors is about 40% larger than the quenched 
value, but that the value of the ratio Ti^/H-i for three light-quark flavors is approximately 
14% lower than the quenched value. According to the velocity-scaling rules, both Ti% and 
H.8 are of order t> 4 (Ref. |1|]). Therefore, we expect M^TisfHi to be of order unity times 
a factor l/(2iV c ), where N c = 3 is the number of QCD colors, to account for the spin and 
color traces in the definition of Tii (Ref. [fTTH ). Our result for M^TisfHi is smaller than 
this estimate by about a factor of three. Our values for Tii and Tis could be used to make 
predictions for the, as yet, unmeasured Xb decay rates. However, as we have mentioned, large 



next-to-leading-order corrections in the perturbation series for those decay rates |L6[ suggest 
that further theoretical progress may be necessary in order to achieve a precise comparison 
with experiment. 

Our results indicate that the quenched approximation yields a poor estimate of the 
NRQCD matrix elements. As we have mentioned, the trends in going from the quenched 
approximation to the physical number of light-quark flavors can be understood in terms of 
a simple picture. The lattice spacing and heavy-quark mass are determined by fitting to 
bottomonium spectroscopy, which probes the wave functions at distances of order l/(Mbv). 
On the other hand, the decay matrix elements sample the wave function and its derivatives 
at much shorter distances, of order a. In the absence of the sea of light quark- ant iquark 
pairs, the strong coupling constant runs too fast, becoming too small at the shorter distance. 
This leads to an underestimate of the decay matrix elements, since the values of the wave 
function and its derivatives depend on the strength of the potential at short distances. We 
note that the S'-wave matrix elements Gi and H 8 both increase by about the same fraction 
in going from the quenched approximation to three light-quark flavors, while the P-wave 
matrix element Hi increases by a larger fraction. This may be because P-wave matrix 
elements depend on the derivative of the wave function at the origin, as opposed to the wave 
function at the origin, and, so, are more sensitive to changes in the strength of the potential 
at short distances. 
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In order to improve the lattice estimates of the matrix elements for S- and P-wave 
bottomonium decay, one would first need more precise measurements of the lattice spacing 
a. A more stringent comparison with experiment would also require a more precise deter- 
mination of the bottom-quark mass Mi,, as well as a calculation to two-loop accuracy of 
the perturbative coefficients that relate the lattice matrix elements to the continuum matrix 
elements.^ In the case of T\ this last improvement is essential to obtain a useful prediction 
of the continuum matrix element. Beyond this it would be valuable to use gauge-field con- 
figurations that have been generated with improved actions. Only when this has been done 
could one justify using NRQCD actions that have been improved to higher orders in v and 
a for the extraction of bottomonium decay matrix elements. 
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Alternatively, one could carry out a non-perturbative (lattice) calculation of these coefficients 
or compute the short-distance decay coefficients in lattice perturbation theory, rather than in 
continuum perturbation theory. 
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